Show the code
#Reading in the West Coast shapefile
wcEEZ <- st_read(here("posts"
, "2022-12-03-geospatial"
, "data"
, "wc_regions_clean.shp"))Ruth Enriquez
December 3, 2022
This work comes from a Geospatial Analyst & Remote Sensing homework assignment. We learned that aquaculture in marine environments holds the promise of being a significant contributor to the world’s food resources, offering a more sustainable protein alternative compared to traditional land-based meat farming.
For this work, I will be identifying which Exclusive Economic Zones (EEZ) on the West Coast of the US are best suited to developing marine aquaculture for several species of oysters.
Based on previous research, we know that oysters needs the following conditions for optimal growth:
sea surface temperature: 11-30°C
depth: 0-70 meters below sea level
Objectives:
combining vector/raster data
resampling raster data
masking raster data
map algebra
I will use average annual sea surface temperature (SST) from the years 2008 to 2012 to characterize the average sea surface temperature within the region. The data we are working with was originally generated from NOAA’s 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1.
To characterize the depth of the ocean I will use the General Bathymetric Chart of the Oceans (GEBCO).3
I will be designating maritime boundaries using Exclusive Economic Zones off of the west coast of US from Marineregions.org.
I will begin by loading all of the necessary packages to analyze the data. Then I will load and validate the data to ensure that it is all on the same coordinate reference system.
::: {.cell}
```{.r .cell-code}
#loading in packages
library(tidyverse)
library(dplyr)
library(here)
library(janitor)
library(tmap)
library(tmaptools)
library(terra)
library(sf)
library(stringr)
# setting file path
setwd(here())
```
:::
read in the shapefile for the West Coast EEZ (wc_regions_clean.shp)
read in SST rasters
average_annual_sst_2008.tifaverage_annual_sst_2009.tifaverage_annual_sst_2010.tifaverage_annual_sst_2011.tifaverage_annual_sst_2012.tifcombine SST rasters into a raster stack
read in bathymetry raster (depth.tif)
check that data are in the same coordinate reference system
Next, I need to process the SST and depth data so that they can be combined. In this case the SST and depth data have slightly different resolutions, extents, and positions. I don’t want to change the underlying depth data, so I will need to resample to match the SST data using the nearest neighbor approach.
find the mean SST from 2008-2012
class : SpatRaster
dimensions : 480, 408, 1 (nrow, ncol, nlyr)
resolution : 0.04165905, 0.04165905 (x, y)
extent : -131.9848, -114.9879, 29.99208, 49.98842 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s) : memory
name : mean
min value : 281.3675
max value : 301.1837
convert SST data from Kelvin to Celsius
crop depth raster to match the extent of the SST raster
note: the resolutions of the SST and depth data do not match
check that the depth and SST match in resolution, extent, and coordinate reference system
class : SpatRaster
dimensions : 480, 408, 1 (nrow, ncol, nlyr)
resolution : 0.04165905, 0.04165905 (x, y)
extent : -131.9848, -114.9879, 29.99208, 49.98842 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s) : memory
name : mean
min value : 8.217548
max value : 28.033722
class : SpatRaster
dimensions : 480, 408, 1 (nrow, ncol, nlyr)
resolution : 0.04165905, 0.04165905 (x, y)
extent : -131.9848, -114.9879, 29.99208, 49.98842 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s) : memory
name : depth
min value : -5468
max value : 4218
#check that the resolution is the same
#resolution written: x, y
#stackSSTmeanC
#resolution : 0.04165905, 0.04165905
#oceanDepthRe
#resolution : 0.04165905, 0.04165905
#check that the extent is the same
#extent written in: xmin, xmax, ymin, ymax
#stackSSTmeanC
#extent: -131.9848, -114.9879, 29.99208, 49.98842
#oceanDepthRe
#extent: -131.9848, -114.9879, 29.99208, 49.98842
#check the crs is the same
#stackSSTmeanC
#coord. ref. : lon/lat WGS 84 (EPSG:4326)
#oceanDepthRe
#coord. ref. : lon/lat WGS 84 (EPSG:4326) In order to find suitable locations for marine aquaculture, I will need to find locations that are suitable in terms of both SST and depth.
reclassify SST and depth data into locations that are suitable for oyster
#identify suitable conditions for oysters
#(taken from data given)...
#sea surface temperature: 11-30°C
#depth: 0-70 meters below sea level
#reclassifying temp
temp <- c(-Inf, 11, NA,
11, 30, 1,
30, Inf, NA)
tempMatrix <- matrix(temp, ncol = 3, byrow = TRUE)
SSTreclass <- classify(SSTcelsius, rcl = tempMatrix, include.lowest = TRUE)
#reclassifying depth
depth <- c(-Inf, -70, NA,
-70, 0, 1,
0, Inf, NA)
depthMatrix <- matrix(depth, ncol = 3, byrow = TRUE)
depthReclass <- classify(oceanDepthRe, rcl = depthMatrix, include.lowest = TRUE)find locations that satisfy both SST and depth conditions
I want to determine the total suitable area within each EEZ in order to rank zones by priority. To do so, we need to find the total area of suitable locations within each EEZ.
#setting up cell size
cellSize <- cellSize(suitableLocations, mask = TRUE, unit = "km", transform = TRUE)
#changing EEZ file to raster
wcRaster <- rasterize(wcEEZ, suitableLocations, field ="rgn")
#creating a mask, testing which variable combos work: wcRaster, suitableLocations, EEZcropSuitable
wcMask <- mask(wcRaster, suitableLocations, inverse = FALSE, updatevalue = NA)
#using zonal example from week 5 lab... applying zonal operations to help us find zones
suitableArea <- zonal(cellSize, wcMask, sum, na.rm = TRUE)
#joining together data
suitableEEZ <- merge(wcEEZ, suitableArea, by = 'rgn') |>
mutate(suitable_area = area, percentage = (suitable_area/area_km2 * 100), .before = geometry)Now that I have results, I need to present them!
Time to create the following maps:
tmap mode set to interactive viewing
#Creating a map for total suitable area by region
tm_basemap("Esri.WorldStreetMap")+
tm_shape(suitableEEZ) +
tm_polygons(col = 'area',
palette = 'BrBG',
alpha = 0.5,
border.col = 'black') +
tm_layout( main.title = "Suitable Areas for Oysters by Region: Total",
main.title.position = 'center',
main.title.size = 0.5,
legend.outside = TRUE) +
tm_scale_bar(position = c("left", "bottom"))+
tm_text("rgn", size = 0.54)#Creating a map for percent suitable area by region
tm_basemap("Esri.WorldStreetMap")+
tm_shape(suitableEEZ) +
tm_polygons(col = 'percentage',
palette = 'BrBG',
alpha = 0.5,
border.col = 'black') +
tm_layout( main.title = "Suitable Areas for Oysters by Region: Percentage",
main.title.position = 'center',
main.title.size = 0.5,
legend.outside = TRUE)+
tm_scale_bar(position = c("left", 'bottom')) +
tm_text("rgn", size = 0.54)Now that I’ve worked through the solution for one group of species, I want to update my workflow to work for other species I am interested in. To do this I will create a function that would allow anyone to reproduce my results for other species. My function will be able to do the following:
I am interested at looking for optimal EEZs for Dungness Crab. To get information on the depth and temperature requirements for Dungess Crab I will go to SeaLifeBase.
# Creating the function
suitableMapFunction <- function(seaSurfaceTempLow, seaSurfaceTempHigh, oceanDepthLow, oceanDepthHigh, speciesName ){
wcEEZ <- st_read(here("posts"
, "2022-12-03-geospatial"
, "data"
, "wc_regions_clean.shp"))
filelist <- list.files(here("posts"
, "2022-12-03-geospatial"
, "data"
, "sst"), full.names = TRUE)
orgStackSST <- rast(filelist)
oceanDepth <- rast(here("posts"
, "2022-12-03-geospatial"
, "data"
, "depth.tif"))
stackSST <- project(orgStackSST, y ="epsg:4326")
stackSSTmean <- mean(stackSST)
SSTcelsius <- stackSSTmean - 273.15
oceanDepthCrop <- crop(oceanDepth, SSTcelsius)
oceanDepthRe <- resample(oceanDepthCrop, SSTcelsius, method = 'near')
temp <- c(-Inf, seaSurfaceTempLow, NA,
seaSurfaceTempLow, seaSurfaceTempHigh, 1,
seaSurfaceTempHigh, Inf, NA)
tempMatrix <- matrix(temp, ncol = 3, byrow = TRUE)
SSTreclass <- classify(SSTcelsius, rcl = tempMatrix, include.lowest = TRUE)
depth <- c(-Inf, oceanDepthLow, NA,
oceanDepthLow, oceanDepthHigh, 1,
oceanDepthHigh, Inf, NA)
depthMatrix <- matrix(depth, ncol = 3, byrow = TRUE)
depthReclass <- classify(oceanDepthRe, rcl = depthMatrix, include.lowest = TRUE)
funC <- function(x, y) {
return(x*y)
}
suitableLocations <-lapp(c(SSTreclass,depthReclass), fun = funC)
cellSize <- cellSize(suitableLocations, mask = TRUE, unit = "km", transform = TRUE)
wcRaster <- rasterize(wcEEZ, suitableLocations, field ="rgn")
wcMask <- mask(wcRaster, suitableLocations, inverse = FALSE, updatevalue = NA)
suitableArea <- zonal(cellSize, wcMask, sum, na.rm = TRUE)
suitableEEZ <- merge(wcEEZ, suitableArea, by = 'rgn') |>
mutate(suitable_area = area, percentage = (suitable_area/area_km2 * 100), .before = geometry)
suitableAreaTotalMap <- tm_shape(suitableEEZ) +
tm_polygons(col = 'area',
palette = 'BrBG',
alpha = 0.5,
border.col = 'black') +
tm_layout(main.title = paste("Suitable Areas for", speciesName, "by Region: Total"),
main.title.position = 'center',
title.size = 0.5,
legend.outside = TRUE) +
tm_scale_bar(position = c("left", "bottom"))+
tm_text("rgn", size = 0.54)
suitableAreaPercentMap <- tm_shape(suitableEEZ) +
tm_polygons(col = 'percentage',
palette = 'BrBG',
alpha = 0.5,
border.col = 'black') +
tm_layout(main.title = paste("Suitable Areas for", speciesName, " by Region: Percentage"),
main.title.position = 'center',
title.size = 0.03,
legend.outside = TRUE)+
tm_scale_bar(position = c("left", 'bottom')) +
tm_text("rgn", size = 0.54)
tmap_arrange(suitableAreaTotalMap, suitableAreaPercentMap, widths = c(.25, .75))
}Reading layer `wc_regions_clean' from data source
`C:\Users\ruthe\Documents\ruthe808.github.io\posts\2022-12-03-geospatial\data\wc_regions_clean.shp'
using driver `ESRI Shapefile'
Simple feature collection with 5 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -129.1635 ymin: 30.542 xmax: -117.097 ymax: 49.00031
Geodetic CRS: WGS 84